LC_thin_compare.R

# This script analyzes the results of LC_processVS_hpc.R

Packages <- c("rstan", "coda", "ggmcmc", "bayesplot", "tidyverse", 
              "loo", "purrr")
suppressMessages(invisible(lapply(Packages, library, character.only=TRUE)))
theme_set(theme_bw())

setwd("C:/Users/tms1044/Desktop/lc_mod")

ggs_autocor <- function (D, family = NA, nLags = 20, greek = FALSE) 
{
  if (!is.na(family)) {
    D <- get_family(D, family = family)
  }
  nIter <- attr(D, "nIteration")

  wc.ac <- D %>% dplyr::group_by(Parameter, Chain) %>% 
    dplyr::do(ac(.$value, nLags))
  f <- ggplot(wc.ac, aes(x=Lag, y=Autocorrelation, colour=as.factor(Chain))) + 
    theme(legend.position="none") + geom_hline(yintercept=0, colour="gray70") +
    geom_point(size=0.5) + geom_line(stat="identity", position="identity") + 
    scale_colour_brewer(type="qual", palette=2) + 
    facet_wrap(~Parameter) + ylim(-.25, 1)
  return(f)
}

thin.seq <- seq(4,20,2)
chn.dir <- "out/"
sum.dir <- "thin_summaries/"

# GGMCMC
beta.fls <- list.files(sum.dir, pattern="out_betas_", full.names=TRUE)
beta.fls <- beta.fls[grepl("_Clim.", beta.fls, fixed=TRUE)]
names(beta.fls) <- beta.fls
beta.fls <- beta.fls[c(7:9,1:6)]
map(beta.fls, readRDS) %>% 
  map2(., beta.fls, ~(ggs_autocor(ggs(.x)) + ggtitle(.y))) %>% print

## $`thin_summaries/out_betas_4_Clim.rds`

## 
## $`thin_summaries/out_betas_6_Clim.rds`

## 
## $`thin_summaries/out_betas_8_Clim.rds`

## 
## $`thin_summaries/out_betas_10_Clim.rds`

## 
## $`thin_summaries/out_betas_12_Clim.rds`

## 
## $`thin_summaries/out_betas_14_Clim.rds`

## 
## $`thin_summaries/out_betas_16_Clim.rds`

## 
## $`thin_summaries/out_betas_18_Clim.rds`

## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`thin_summaries/out_betas_20_Clim.rds`

map(beta.fls, readRDS) %>% 
  map2(., beta.fls, ~(ggs_density(ggs(.x)) + ggtitle(.y) +
                        facet_wrap(~Parameter, scales="free"))) %>% print

## $`thin_summaries/out_betas_4_Clim.rds`

## 
## $`thin_summaries/out_betas_6_Clim.rds`

## 
## $`thin_summaries/out_betas_8_Clim.rds`

## 
## $`thin_summaries/out_betas_10_Clim.rds`

## 
## $`thin_summaries/out_betas_12_Clim.rds`

## 
## $`thin_summaries/out_betas_14_Clim.rds`

## 
## $`thin_summaries/out_betas_16_Clim.rds`

## 
## $`thin_summaries/out_betas_18_Clim.rds`

## 
## $`thin_summaries/out_betas_20_Clim.rds`

map(beta.fls, readRDS) %>% 
  map2(., beta.fls, ~(ggs_traceplot(ggs(.x)) + ggtitle(.y) + 
                        facet_wrap(~Parameter, scales="free"))) %>% print

## $`thin_summaries/out_betas_4_Clim.rds`

## 
## $`thin_summaries/out_betas_6_Clim.rds`

## 
## $`thin_summaries/out_betas_8_Clim.rds`

## 
## $`thin_summaries/out_betas_10_Clim.rds`

## 
## $`thin_summaries/out_betas_12_Clim.rds`

## 
## $`thin_summaries/out_betas_14_Clim.rds`

## 
## $`thin_summaries/out_betas_16_Clim.rds`

## 
## $`thin_summaries/out_betas_18_Clim.rds`

## 
## $`thin_summaries/out_betas_20_Clim.rds`

map(beta.fls, readRDS) %>% 
  map2(., beta.fls, ~(ggs_caterpillar(ggs(.x)) + ggtitle(.y))) %>% print

## $`thin_summaries/out_betas_4_Clim.rds`

## 
## $`thin_summaries/out_betas_6_Clim.rds`

## 
## $`thin_summaries/out_betas_8_Clim.rds`

## 
## $`thin_summaries/out_betas_10_Clim.rds`

## 
## $`thin_summaries/out_betas_12_Clim.rds`

## 
## $`thin_summaries/out_betas_14_Clim.rds`

## 
## $`thin_summaries/out_betas_16_Clim.rds`

## 
## $`thin_summaries/out_betas_18_Clim.rds`

## 
## $`thin_summaries/out_betas_20_Clim.rds`

map(beta.fls, readRDS) %>% 
  map2(., beta.fls, ~(ggs_geweke(ggs(.x)) + ggtitle(.y))) %>% print

## $`thin_summaries/out_betas_4_Clim.rds`

## 
## $`thin_summaries/out_betas_6_Clim.rds`

## 
## $`thin_summaries/out_betas_8_Clim.rds`

## 
## $`thin_summaries/out_betas_10_Clim.rds`

## 
## $`thin_summaries/out_betas_12_Clim.rds`

## 
## $`thin_summaries/out_betas_14_Clim.rds`

## 
## $`thin_summaries/out_betas_16_Clim.rds`

## 
## $`thin_summaries/out_betas_18_Clim.rds`

## 
## $`thin_summaries/out_betas_20_Clim.rds`

map(beta.fls, readRDS) %>% 
  map2(., beta.fls, ~(ggs_Rhat(ggs(.x)) + ggtitle(.y))) %>% print

## $`thin_summaries/out_betas_4_Clim.rds`

## 
## $`thin_summaries/out_betas_6_Clim.rds`

## 
## $`thin_summaries/out_betas_8_Clim.rds`

## 
## $`thin_summaries/out_betas_10_Clim.rds`

## 
## $`thin_summaries/out_betas_12_Clim.rds`

## 
## $`thin_summaries/out_betas_14_Clim.rds`

## 
## $`thin_summaries/out_betas_16_Clim.rds`

## 
## $`thin_summaries/out_betas_18_Clim.rds`

## 
## $`thin_summaries/out_betas_20_Clim.rds`

# effective sample size
par(mfrow=c(3,3))
map(beta.fls, readRDS) %>% 
  walk2(., thin.seq, ~hist(effectiveSize(.x), main=.y, xlab="n_eff"))

par(mfrow=c(3,3))
map(beta.fls, readRDS) %>% 
  walk2(., thin.seq, ~hist(effectiveSize(.x)/(5000*6/.y), main=.y,
                           xlab="n_eff/n_iter", xlim=c(0,1.5)))

par(mfrow=c(1,1))